home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.4 / XPM / BDY.C next >
C/C++ Source or Header  |  1999-09-11  |  16KB  |  602 lines
  1. /* 
  2.  * bdy.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * BDY.C
  11.  *
  12.  * functions implementating C.T.Zahn's algorithm for boundary
  13.  * detection and representation
  14.  *
  15.  */
  16. #include <stdio.h>
  17. #include <malloc.h>
  18. #include <stdlib.h>
  19. #include <math.h>
  20. #include "ip.h"
  21. #include "bdy.h"
  22.  
  23. #define    WINDOW_ROWS    3
  24. #define NDIRECTIONS    8
  25.  
  26. #define    THRESH        1
  27. #define MIN_POLY_SIZE    1
  28. #define    ZERO        0
  29. #define ROOT2        1.414213562
  30.  
  31. #define    ON        1
  32. #define    OFF        0
  33. #undef    SHOW_CURV_PT
  34.  
  35. #define    DISPL_PAGE    1
  36. #define    RESET        ON
  37.  
  38. #define IN_LIST        1
  39. #define    OUT_LIST    0
  40. #define SHOW_CURV_PT
  41.  
  42.  
  43. /* global variables */
  44. extern struct curv_points *curv_head_in[NDIRECTIONS];
  45. extern struct curv_points *curv_head_out[NDIRECTIONS];
  46. extern struct curv_points *curv_tail_in[NDIRECTIONS];
  47. extern struct curv_points *curv_tail_out[NDIRECTIONS];
  48.  
  49. extern struct polygon *poly_head_ptr;
  50. extern struct polygon *poly_tail_ptr;
  51.  
  52. extern float *delta_phik, *delta_lk;
  53. extern long ncp;
  54.  
  55. int page = DISPL_PAGE;
  56.  
  57. unsigned int d1, d2, d3, d4, d5, d6, b1, b2, b3, b4, b5, b6;
  58. unsigned nbytes = 1;
  59.  
  60.  
  61. /*
  62.  * generate curvature points using Zahn algorithm
  63.  * (ref: C. T. Zahn, SLAC report No. 70, Dec. 1966)
  64.  */
  65. void
  66. curvature_points (unsigned char *window[], int ir, int jmax, int imax)
  67. {
  68.   int jc;
  69.  
  70.   for (jc = 0; jc < jmax - 2; jc++) {
  71.     d1 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  72.     d2 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  73.     d3 = *(*(window + (ir % WINDOW_ROWS)) + jc + 2);
  74.     d4 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  75.     d5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  76.     d6 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 2);
  77.     horizontal (ir, jc);
  78.   }
  79.   for (jc = 0; jc < jmax - 1; jc++) {
  80.     b1 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  81.     b2 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  82.     b3 = *(*(window + ((ir + 2) % WINDOW_ROWS)) + jc + 1);
  83.     b4 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  84.     b5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  85.     b6 = *(*(window + ((ir + 2) % WINDOW_ROWS)) + jc);
  86.     vertical (ir, jc);
  87.   }
  88.   if (ir == imax - WINDOW_ROWS) {
  89.     ir++;
  90.     for (jc = 0; jc < jmax - 2; jc++) {
  91.       d1 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  92.       d2 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  93.       d3 = *(*(window + (ir % WINDOW_ROWS)) + jc + 2);
  94.       d4 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  95.       d5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  96.       d6 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 2);
  97.       horizontal (ir, jc);
  98.     }
  99.   }
  100. }
  101.  
  102.  
  103. /*
  104.  * check 3x2 horizontal array (in a binary image) and determine 
  105.  * direction codes for incoming and outgoing contour segment
  106.  */
  107. int
  108. horizontal (i, j)
  109.      int i, j;
  110. {
  111.   int in, out;
  112.   int err_flag = -1;
  113.  
  114.   if (d2 >= THRESH) {
  115.     if (d5 == ZERO) {
  116.       if (d6 >= THRESH)
  117.         in = 3;
  118.       else if (d3 >= THRESH)
  119.         in = 4;
  120.       else
  121.         in = 5;
  122.  
  123.       if (d4 >= THRESH)
  124.         out = 5;
  125.       else if (d1 >= THRESH)
  126.         out = 4;
  127.       else
  128.         out = 3;
  129.     }
  130.     else
  131.       return (err_flag);
  132.   }
  133.   else {
  134.     if (d5 >= THRESH) {
  135.       if (d1 >= THRESH)
  136.         in = 7;
  137.       else if (d4 >= THRESH)
  138.         in = 0;
  139.       else
  140.         in = 1;
  141.  
  142.       if (d3 >= THRESH)
  143.         out = 1;
  144.       else if (d6 >= THRESH)
  145.         out = 0;
  146.       else
  147.         out = 7;
  148.     }
  149.     else
  150.       return (err_flag);
  151.   }
  152.  
  153. /*
  154.  * curvature point found! 
  155.  * allocate memory and store it in linked list curv_points structure 
  156.  */
  157.   if (in != out) {
  158.  
  159. #ifdef SHOW_CURV_PT
  160.     printf ("   horizontal:  ir= %d, jc= %d, in= %d, out= %d\n",
  161.             2 * (i + 1), 2 * (j + 1) + 1, in, out);
  162. #endif
  163.     if (curv_head_in[in] == NULL) {
  164.       curv_head_in[in] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  165.       curv_head_in[in]->prev = NULL;
  166.       curv_tail_in[in] = curv_head_in[in];
  167.     }
  168.     else {
  169.       curv_tail_in[in]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  170.       curv_tail_in[in]->next->prev = curv_tail_in[in];
  171.       curv_tail_in[in] = curv_tail_in[in]->next;
  172.     }
  173.     if (curv_head_out[out] == NULL) {
  174.       curv_head_out[out] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  175.       curv_head_out[out]->prev = NULL;
  176.       curv_tail_out[out] = curv_head_out[out];
  177.     }
  178.     else {
  179.       curv_tail_out[out]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  180.       curv_tail_out[out]->next->prev = curv_tail_out[out];
  181.       curv_tail_out[out] = curv_tail_out[out]->next;
  182.     }
  183.     ncp++;
  184.  
  185.     curv_tail_in[in]->same_point = curv_tail_out[out];
  186.     curv_tail_in[in]->next = NULL;
  187.     curv_tail_out[out]->next = NULL;
  188.     curv_tail_in[in]->xn = curv_tail_out[out]->xn = 2 * (j + 1) + 1;
  189.     curv_tail_in[in]->yn = curv_tail_out[out]->yn = 2 * (i + 1);
  190.     curv_tail_in[in]->edge_in = curv_tail_out[out]->edge_in = in;
  191.     curv_tail_in[in]->edge_out = curv_tail_out[out]->edge_out = out;
  192.   }
  193.   return (1);
  194. }
  195.  
  196.  
  197. /*
  198.  * check 2x3 vertical array (in a binary image) and determine 
  199.  * direction codes for incoming and outgoing contour segment
  200.  */
  201. int
  202. vertical (i, j)
  203.      int i, j;
  204. {
  205.   int in, out;
  206.   int err_flag = -1;
  207.  
  208.   if (b2 >= THRESH) {
  209.     if (b5 == ZERO) {
  210.       if (b6 >= THRESH)
  211.         in = 1;
  212.       else if (b3 >= THRESH)
  213.         in = 2;
  214.       else
  215.         in = 3;
  216.  
  217.       if (b4 >= THRESH)
  218.         out = 3;
  219.       else if (b1 >= THRESH)
  220.         out = 2;
  221.       else
  222.         out = 1;
  223.     }
  224.     else
  225.       return (err_flag);
  226.   }
  227.   else {
  228.     if (b5 >= THRESH) {
  229.       if (b1 >= THRESH)
  230.         in = 5;
  231.       else if (b4 >= THRESH)
  232.         in = 6;
  233.       else
  234.         in = 7;
  235.  
  236.       if (b3 >= THRESH)
  237.         out = 7;
  238.       else if (b6 >= THRESH)
  239.         out = 6;
  240.       else
  241.         out = 5;
  242.     }
  243.     else
  244.       return (err_flag);
  245.   }
  246.  
  247. /*
  248.  * curvature point found! 
  249.  * allocate memory and store it in linked list curv_points structure 
  250.  */
  251.   if (in != out) {
  252.  
  253. #ifdef SHOW_CURV_PT
  254.     printf ("  vertical:  ir= %d, jc= %d, in= %d, out= %d\n",
  255.             2 * (i + 1) + 1, 2 * (j + 1), in, out);
  256. #endif
  257.  
  258.     if (curv_head_in[in] == NULL) {
  259.       curv_head_in[in] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  260.       curv_head_in[in]->prev = NULL;
  261.       curv_tail_in[in] = curv_head_in[in];
  262.     }
  263.     else {
  264.       curv_tail_in[in]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  265.       curv_tail_in[in]->next->prev = curv_tail_in[in];
  266.       curv_tail_in[in] = curv_tail_in[in]->next;
  267.     }
  268.     if (curv_head_out[out] == NULL) {
  269.       curv_head_out[out] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  270.       curv_head_out[out]->prev = NULL;
  271.       curv_tail_out[out] = curv_head_out[out];
  272.     }
  273.     else {
  274.       curv_tail_out[out]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  275.       curv_tail_out[out]->next->prev = curv_tail_out[out];
  276.       curv_tail_out[out] = curv_tail_out[out]->next;
  277.     }
  278.     ncp++;
  279.  
  280.     curv_tail_in[in]->same_point = curv_tail_out[out];
  281.     curv_tail_in[in]->next = NULL;
  282.     curv_tail_out[out]->next = NULL;
  283.     curv_tail_in[in]->xn = curv_tail_out[out]->xn = 2 * (j + 1);
  284.     curv_tail_in[in]->yn = curv_tail_out[out]->yn = 2 * (i + 1) + 1;
  285.     curv_tail_in[in]->edge_in = curv_tail_out[out]->edge_in = in;
  286.     curv_tail_in[in]->edge_out = curv_tail_out[out]->edge_out = out;
  287.   }
  288.   return (1);
  289. }
  290.  
  291.  
  292. /*
  293.  * link up all curvature points into closed polygons
  294.  */
  295. void
  296. linkage (Image * imgIn, Image * imgOut, int value)
  297. {
  298.   struct polygon *dummy_poly_ptr;
  299.   struct curv_points *first_point;
  300.   float *d_phi_head, *d_l_head;
  301.   int m, c;
  302.   char in_buf[IN_BUF_LEN];
  303.  
  304.   do {
  305.     d_phi_head = delta_phik;
  306.     d_l_head = delta_lk;
  307.     for (m = 0; m < NDIRECTIONS; m++)
  308.       if ((first_point = curv_head_out[m]) != NULL) {
  309.         poly (first_point);
  310.         break;
  311.       }
  312.   } while (first_point != NULL);
  313.  
  314.   if (ncp != 0)
  315.     printf ("...something wrong, didn't find all polygons!\n");
  316.   if (poly_head_ptr == NULL) {
  317.     printf ("...no objects found to analyze!\n");
  318.     return;
  319.   }
  320.  
  321. /*
  322.  * polygon analysis
  323.  */
  324.   printf ("...analyze a boundary?(y/n) -- ");
  325.   while ((c = readlin (in_buf)) == 'y') {
  326.     dummy_poly_ptr = select_poly (poly_head_ptr, imgOut, value);
  327.     printf ("\n...analyze another boundary?(y/n)");
  328.   }
  329. }
  330.  
  331. /*
  332.  * create a closed polygon given a starting curvature point 
  333.  */
  334. void
  335. poly (first_point)
  336.      struct curv_points *first_point;
  337. {
  338.   struct curv_points *next_point, *prev_point, *match_in_list ();
  339.  
  340.   if (poly_head_ptr == NULL) {
  341.     poly_head_ptr = (struct polygon *) calloc (nbytes, sizeof (struct polygon));
  342.     poly_tail_ptr = poly_head_ptr;
  343.     poly_head_ptr->prev_poly = NULL;
  344.   }
  345.   else {
  346.     poly_tail_ptr->next_poly = (struct polygon *) calloc (nbytes, sizeof (struct polygon));
  347.     poly_tail_ptr->next_poly->prev_poly = poly_tail_ptr;
  348.     poly_tail_ptr = poly_tail_ptr->next_poly;
  349.   }
  350.   poly_tail_ptr->next_poly = NULL;
  351.  
  352.   prev_point = first_point;
  353.   next_point = match_in_list (prev_point);
  354.   poly_tail_ptr->d_phi_ptr = delta_phik;
  355.   poly_tail_ptr->d_l_ptr = delta_lk;
  356.   poly_tail_ptr->poly_points = 0;
  357.   poly_tail_ptr->total_delta_phi = (float) 0;
  358.   poly_tail_ptr->first_x = first_point->xn;
  359.   poly_tail_ptr->first_y = first_point->yn;
  360.   poly_tail_ptr->first_edge_out = first_point->edge_out;
  361.   create_edge (prev_point, next_point);
  362.   delete_list (prev_point, OUT_LIST);
  363.   prev_point = next_point->same_point;
  364.   delete_list (next_point, IN_LIST);
  365.   ncp--;
  366.   poly_tail_ptr->poly_points++;
  367.   poly_tail_ptr->total_delta_phi = poly_tail_ptr->total_delta_phi + *delta_phik;
  368.   delta_lk++;
  369.   delta_phik++;
  370.   do {
  371.     next_point = match_in_list (prev_point);
  372.     create_edge (prev_point, next_point);
  373.     delete_list (prev_point, OUT_LIST);
  374.     prev_point = next_point->same_point;
  375.     delete_list (next_point, IN_LIST);
  376.     ncp--;
  377.     poly_tail_ptr->poly_points++;
  378.     poly_tail_ptr->total_delta_phi = poly_tail_ptr->total_delta_phi
  379.       + *delta_phik;
  380.     delta_lk++;
  381.     delta_phik++;
  382.  
  383.     if (ncp < 0) {
  384.       printf ("Something wrong, points<0\n");
  385.       exit (1);
  386.     }
  387.   }
  388.   while (prev_point != first_point);
  389.  
  390.   if (poly_tail_ptr->poly_points <= MIN_POLY_SIZE) {
  391.     if (poly_tail_ptr->prev_poly == NULL)
  392.       poly_head_ptr = NULL;
  393.     else {
  394.       poly_tail_ptr = poly_tail_ptr->prev_poly;
  395.       free (poly_tail_ptr->next_poly);
  396.       poly_tail_ptr->next_poly = NULL;
  397.     }
  398.   }
  399. }
  400.  
  401. /*
  402.  * given a point in the out-list, find the next point
  403.  * in the in-list with edge-out(out-list) = edge-in(in-list)
  404.  */
  405. struct curv_points *
  406. match_in_list (current_point)
  407.      struct curv_points *current_point;
  408. {
  409.   struct curv_points *tail_point, *head_point;
  410.   int direction;
  411.   int xcurr, ycurr, xnext, ynext;
  412.  
  413.   xcurr = current_point->xn;
  414.   ycurr = current_point->yn;
  415.   printf ("match_in_list: current point = (%d, %d)\n", xcurr, ycurr);
  416.   direction = current_point->edge_out;
  417.   tail_point = curv_tail_in[direction];
  418.   head_point = curv_head_in[direction];
  419.   while ((head_point != NULL) || (tail_point != NULL)) {
  420.     switch (direction) {
  421.     case 0:
  422.       ynext = head_point->yn;
  423.       xnext = head_point->xn;
  424.       if ((ycurr == ynext) && (xnext > xcurr))
  425.         return (head_point);
  426.       break;
  427.     case 4:
  428.       ynext = tail_point->yn;
  429.       xnext = tail_point->xn;
  430.       if ((ycurr == ynext) && (xnext < xcurr))
  431.         return (tail_point);
  432.       break;
  433.     case 2:
  434.       ynext = tail_point->yn;
  435.       xnext = tail_point->xn;
  436.       if ((xcurr == xnext) && (ynext < ycurr))
  437.         return (tail_point);
  438.       break;
  439.     case 6:
  440.       ynext = head_point->yn;
  441.       xnext = head_point->xn;
  442.       if ((xcurr == xnext) && (ynext > ycurr))
  443.         return (head_point);
  444.       break;
  445.     case 1:
  446.       ynext = tail_point->yn;
  447.       xnext = tail_point->xn;
  448.       if (((xcurr + ycurr) == (xnext + ynext)) &&
  449.           (xnext > xcurr) && (ynext < ycurr))
  450.         return (tail_point);
  451.       break;
  452.     case 5:
  453.       ynext = head_point->yn;
  454.       xnext = head_point->xn;
  455.       if (((xcurr + ycurr) == (xnext + ynext)) &&
  456.           (xnext < xcurr) && (ynext > ycurr))
  457.         return (head_point);
  458.       break;
  459.     case 3:
  460.       ynext = tail_point->yn;
  461.       xnext = tail_point->xn;
  462.       if (((xcurr - ycurr) == (xnext - ynext)) &&
  463.           (xnext < xcurr) && (ynext < ycurr))
  464.         return (tail_point);
  465.       break;
  466.     case 7:
  467.       ynext = head_point->yn;
  468.       xnext = head_point->xn;
  469.       if (((xcurr - ycurr) == (xnext - ynext)) &&
  470.           (xnext > xcurr) && (ynext > ycurr))
  471.         return (head_point);
  472.       break;
  473.     }
  474.     tail_point = tail_point->prev;
  475.     head_point = head_point->next;
  476.   }
  477.   printf ("...no match found in outlist[%d]\n", direction);
  478.   exit (1);
  479. }
  480.  
  481. /*
  482.  * delete curvature point from either in-list or out-list
  483.  */
  484. void
  485. delete_list (current_point, list)
  486.      struct curv_points *current_point;
  487.      int list;
  488. {
  489.   int direction;
  490.  
  491.   if (current_point == NULL)
  492.     printf ("...current point can't be NULL!\n");
  493.  
  494.   if ((current_point->prev == NULL) && (current_point->next == NULL)) {
  495.     if (list == IN_LIST) {
  496.       direction = current_point->edge_in;
  497.       curv_head_in[direction] = NULL;
  498.       curv_tail_in[direction] = NULL;
  499.     }
  500.     else {
  501.       direction = current_point->edge_out;
  502.       curv_head_out[direction] = NULL;
  503.       curv_tail_out[direction] = NULL;
  504.     }
  505.   }
  506.   else if (current_point->prev == NULL) {
  507.     if (list == IN_LIST) {
  508.       direction = current_point->edge_in;
  509.       curv_head_in[direction] = current_point->next;
  510.     }
  511.     else {
  512.       direction = current_point->edge_out;
  513.       curv_head_out[direction] = current_point->next;
  514.     }
  515.     current_point->next->prev = NULL;
  516.   }
  517.   else if (current_point->next == NULL) {
  518.     if (list == IN_LIST) {
  519.       direction = current_point->edge_in;
  520.       curv_tail_in[direction] = current_point->prev;
  521.     }
  522.     else {
  523.       direction = current_point->edge_out;
  524.       curv_tail_out[direction] = current_point->prev;
  525.     }
  526.     current_point->prev->next = NULL;
  527.   }
  528.   else {
  529.     current_point->prev->next = current_point->next;
  530.     current_point->next->prev = current_point->prev;
  531.   }
  532.   free (current_point);
  533. }
  534.  
  535.  
  536. void
  537. print_curv ()
  538. {
  539.   int x, y;
  540.   struct curv_points *temp_head;
  541.  
  542.   for (x = 0; x < NDIRECTIONS; x++) {
  543.     y = 0;
  544.     temp_head = curv_head_in[x];
  545.     while (temp_head != NULL) {
  546.       y++;
  547.       temp_head = temp_head->next;
  548.     }
  549.     if (y != 0)
  550.       printf ("inlist[%d] = %d\n", x, y);
  551.   }
  552.  
  553.   for (x = 0; x < NDIRECTIONS; x++) {
  554.     y = 0;
  555.     temp_head = curv_head_out[x];
  556.     while (temp_head != NULL) {
  557.       y++;
  558.       temp_head = temp_head->next;
  559.     }
  560.     if (y != 0)
  561.       printf ("outlist[%d] = %d\n", x, y);
  562.   }
  563. }
  564.  
  565.  
  566. /*
  567.  * create an edge from old-point to new-point
  568.  */
  569. void
  570. create_edge (old_point, new_point)
  571.      struct curv_points *old_point, *new_point;
  572. {
  573.   int phi;
  574.  
  575.   switch (new_point->edge_in) {
  576.   case 0:
  577.   case 4:
  578.     *delta_lk = (float) abs (new_point->xn - old_point->xn);
  579.     phi = new_point->edge_out - new_point->edge_in;
  580.     if (phi == 7)
  581.       phi = -1;
  582.     *delta_phik = (float) phi;
  583.     break;
  584.   case 2:
  585.   case 6:
  586.     *delta_lk = (float) abs (new_point->yn - old_point->yn);
  587.     *delta_phik = (float) (new_point->edge_out - new_point->edge_in);
  588.     break;
  589.   default:
  590.     *delta_lk = (float) abs (new_point->xn - old_point->xn) * (float) ROOT2;
  591.     phi = new_point->edge_out - new_point->edge_in;
  592.     if (phi == -7)
  593.       phi = 1;
  594.     else if (phi == -6)
  595.       phi = 2;
  596.     else if (phi == 6)
  597.       phi = -2;
  598.     *delta_phik = (float) phi;
  599.     break;
  600.   }
  601. }
  602.